An alternative to the plasma emission model: Particle-In-Cell, self-consistent 
electromagnetic wave emission simulations of solar type III radio bursts 
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High-resolution (sub-Debye length grid size and 10000 particle species per cell), 1.5D Particle-in- 
Cell, relativistic, fully electromagnetic simulations are used to model electromagnetic wave emission 
generation in the context of solar type III radio bursts. The model studies generation of electromag- 
netic waves by a super-thermal, hot beam of electrons injected into a plasma thread that contains 
uniform longitudinal magnetic field and a parabolic density gradient. In effect, a single magnetic 
line connecting Sun to earth is considered, for which five cases are studied, (i) We find that the 
physical system without a beam is stable and only low amplitude level electromagnetic drift waves 
(noise) are excited, (ii) The beam injection direction is controlled by setting either longitudinal 
or oblique electron initial drift speed, i.e. by setting the beam pitch angle (the angle between the 
beam velocity vector and the direction of background magnetic field). In the case of zero pitch 
angle i.e. when vt • E± = 0, the beam excites only electrostatic, standing waves, oscillating at local 
plasma frequency, in the beam injection spatial location, and only low level electromagnetic drift 
wave noise is also generated, (iii) In the case of oblique beam pitch angles, i.e. when vt ■ E± 0, 
again electrostatic waves with same properties are excited. However, now the beam also generates 
the electromagnetic waves with the properties commensurate to type III radio bursts. The latter is 
evidenced by the wavelet analysis of transverse electric field component, which shows that as the 
beam moves to the regions of lower density and hence lower plasma frequency, frequency of the elec- 
tromagnetic waves drops accordingly, (iv) When the density gradient is removed, an electron beam 
with an oblique pitch angle still generates the electromagnetic radiation. However, in the latter case 
no frequency decrease is seen, (v) Since in most of the presented results the ratio of electron plasma 
and cyclotron frequencies is close to unity near the beam injection location, in order to prove that 
the generated by the non-zero pitch angle beam electromagnetic emission oscillates at the plasma 
frequency, we also consider a case when the magnetic field (and the cyclotron frequency) is ten 
times smaller. Within the limitations of the model, the study presents the first attempt to produce 
synthetic (simulated) dynamical spectrum of the type III radio bursts in the fully kinetic plasma 
model. The latter is based on 1.5D non-zero pitch angle (non-gyrotropic) electron beam, that is an 
alternative to the plasma emission classical mechanism for which two spatial dimensions are needed. 

PACS numbers: 52.40. Mj;52.25.Os;96.60.Tf;52.35.Qz;52.65.Rr 



I. INTRODUCTION 



The type III solar radio bursts are known to be gener- 
ated by the super-thermal beams of electrons that travel 
away from the Sun on open magnetic field lines [U-Q- 
The beams are likely to be manifestations of magnetic 
reconnection which, in turn, is driven by solar flares. 
However, flares can also drive dispersive Alfven waves 
which also can serve as a source of super-thermal beams. 
In this work we do not focus on a question what is an 
actual source of a beam. Instead, we consider a situation 
when a hot 6 x 10 6 K, super-thermal (vb = 0.5c) beam is 
injected into a cool 3 x 10 5 K, Maxwellian plasma with 
parabolically decreasing density gradient, along an open 
magnetic field line with B = 30 gauss. The latter mimics 
a magnetic field line that connects Sun to Earth. There 
is large body of work done from the observational, mod- 
elling and theoretical viewpoint. We refer the interested 
reader to appropriate reviews 043 and to references in 
Ref. Also Introduction section of Ref. [§] provides 
a good, critical overview of possible mechanisms which 
generate the type III burst electromagnetic (EM) radia- 



tion. In brief, there are three categories of models of type 
III solar radio bursts: (i) Quasilinear theory that uses 
kinetic Fokker-Planck type equation for describing the 
dynamics of an electron beam, in conjunction with the 
spectral energy density evolutionary equations for Lang- 
muir and ion-sound waves. In these models the spectral 
energy density of the Langmuir wave packets (that are ex- 
cited by the bump-on-tail unstable beam) travels along 
the open magnetic field lines with a constant speed and 
this is despite the quasilinear relaxation (formation of a 
plateau in the electron distribution function). This im- 
plies some sort of beam marginal stabilisation (Tol - tl5| . 
Some models also include EM emission into the quasilin- 
ear theory based on so-called drift approximation [H|jl8| , 
where nonlinear beam stabilisation during its propaga- 
tion (so called free streaming) is based on Langmuir-ion 
acoustic wave coupling via ponder-motive force and EM 
emission is prescribed by a power law of the beam to am- 
bient plasma number density ratio. Such models can be 
used to construct and constrain the observed dynamical 
spectra physical parameters, (ii) Stochastic growth the- 
ory p^. l20j|. where density irregularities produce a ran- 
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dom growth, in such a way that Langmuir waves are gen- 
erated stochastically and quasilinear interactions within 
the Langmuir clumps cause the beam to fluctuate about 
marginal stability. Such models can be also used for di- 
rect comparison with the solar type III bursts [2lJ. (iii) 
Direct, kinetic simulation approach of type III bursts 
[23 - [25j | to this date used Particle-In-Cell (PIC) numerical 
method. These models mainly focus on the understand- 
ing of basic physics rather than direct comparison with 
the observations. This is due to the size of simulation do- 
main of the models being too small (only few 1000 Debye 
lengths which is roughly l/10 10 th of 1 AU). 

In Ref.jH we have used 1.5D Vlasov-Maxwell simula- 
tions to model EM emission generation in a fully self- 
consistent plasma kinetic model in the solar physics con- 
text. The simulations presented the generation of EM 
emission by the beam-generated Langmuir waves and 
Larmor drift instability in a plasma thread that connects 
the Sun to Earth with the spatial scales compressed ap- 
propriately. We investigated the effects of spatial den- 
sity gradients on the generation of EM radiation. In the 
case without an electron beam, we found that the inho- 
mogeneous plasma with a uniform background magnetic 
field directed transverse to the density gradient is ape- 
riodically unstable to the Larmor-drift instability. The 
latter produced a novel effect of generation of EM emis- 
sion at plasma frequency. The main results of Ref. 8] 
can be summarised as following: In the case without an 
electron beam, the induced perturbations consist of two 
parts: (i) non-escaping Langmuir-type oscillations which 
are localised in the regions of density inhomogeneity, and 
are highly filamentary, with the period of appearance of 
the filaments close to electron plasma frequency in the 
dense regions; and (ii) escaping EM radiation with phase 
speeds close to the speed of light. When we removed the 
density gradient (i.e. which then makes the plasma sta- 
ble to Larmor-drift instability) and a low density, super- 
thermal, hot beam is injected along the domain, in the 
direction perpendicular to the magnetic field (as in so- 
lar coronal magnetic traps which tend to accelerate the 
particles in the direction perpendicular to the magnetic 
field [26j|), the electron beam quasilinear relaxation gen- 
erates non-escaping Langmuir type oscillations which in 
turn generate escaping EM radiation. We found that 
in the spatial location where the beam is injected, the 
standing waves, oscillating at the plasma frequency, are 
excited. It was suggested that these can be used to in- 
terpret the horizontal strips (the narrow-band line emis- 
sion) observed in some dynamical spectra (27[. We have 
also corroborated quasilinear theory predictions: (i) the 
electron free streaming and (ii) the beam long relaxation 
time, in accord with the analytic expressions. We also 
studied the interplay of Larmor-drift instability and the 
generation of EM emission by the Langmuir waves by 
considering dense electron beam in the Larmor-drift un- 
stable (inhomogeneous) plasma. This enabled us to study 
the deviations from the quasilinear theory. 

In the present study we consider a situation that is 



more relevant to type III radio bursts. The VALIS, 1.5D 
Vlasov-Maxwell code used in Ref.[8j did not allow us to 
set the background magnetic field along the physical do- 
main (along rc-axis) because it only solves for (E x ,E y , 0) 
and (0,0, B z ) EM field components. For this reason the 
results of Ref. [|| were affected by the Larmor-drift insta- 
bility. Thus, they were more applicable to interpreting 
the narrow-band line emission. Because, we considered 
spatially ID situation in Ref.[8j we had to set only one 
grid in the ignorable y-direction. In the latter case, in 
the VALIS code a fluid equation is used to update the 
fluid velocity in the y- direction rather than Vlasov's 
equation. Since there is no pressure gradient in the y- 
direction, which is ignorable, the temperature plays no 
role. Thus, it was not possible to set the electron beam 
velocity/momentum y-component. In turn, the only op- 
tion to set finite k± for the beam (which is a requirement 
to excite EM wave, i.e. to couple the beam to EM emis- 
sion) was to set small, finite perpendicular background 
magnetic field B z q. In effect, having finite k± for the 
beam means that v\, ■ Ej_ ^= 0, i.e. the beam veloc- 
ity/momentum vector has a projection on the transverse 
EM component. Only in this case (in the case of non- 
zero pitch angle) the beam can couple to EM wave [28l. 
Another way to look at the case considered in Ref. 8] 
is to realise that vj, x B z q = E (because the resistivity 
is zero and plasma beta is small, electrons tend to be 
magnetised, 'frozen-into' plasma). In this case the lat- 
ter vector product gives v x B z q = E±_ = E y . This is of 
course plausible in the solar corona, as the magnetic field 
has all three components, but then the situation does 
not adequately describe type III bursts in which the elec- 
tron beams are believed to propagate along (not across!) 
the magnetic field lines. For these reasons, now using 
EPOCH Particle-In-Cell (PIC) code, which can update 
all EM components, and allows to set non-zero Bq x , we 
can (a) suppress the Larmor drift instability; (b) set a fi- 
nite electron beam velocity y-component at t = (hence 
to have Vb ■ E± ^ 0) which can readily excite EM waves; 
(c) hence consider the physical system that adequately 
describes type III radio burst magnetic field geometry. 

We would like to stress that the EM emission in our 
model is different from the classical plasma emission 
mechanism. We elaborate on the difference in the con- 
clusions section. 

The paper is organised as following: In sections Ill.a- 
III. c and IILe we consider the inhomogeneous density 
plasma, with the density profile commensurate to type 
III radio bursts. In Section III. a we present the results 
of an equilibrium test run where the initial conditions 
described below are evolved without imposing an elec- 
tron beam; In Section Ill.b we inject an electron beam 
strictly along the background magnetic field, Bq x , with 
Po x = 0.5jm e c (7 = 1/vi — 0.5 2 = 1.155 everywhere) 
and po y = 0, thus setting Vb ■ E± = 0, in turn expecting 
that only electrostatic (ES) plasma waves to be excited. 
In Section III.c the electron beam is injected at an oblique 
angle with pq x = 0.5"/m e c and po y = 0.5"/m e c, thus set- 
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ting Vb ■ E± 0, in turn expecting that EM waves to 
be excited. For Section III.c we produce the synthetic 
(simulated) dynamical spectrum for the EM waves by 
studying the behaviour of frequency of the EM emission 
generated by the beam as a function of time, which is ex- 
pected to decrease, as the beam movies into the regions 
of decreased density (hence decreased plasma frequency 
u> pe ). Because as we will show below in our model the 
EM emission has frequency close to the plasma frequency, 
u!. pe , the fact that u) pe oc y/nl ensures that the frequency 
of the EM emission decreases in time as the non-zero 
pitch angle electron beam moves towards the regions of 
progressively smaller background plasma density. In Sec- 
tion IILd we consider the situation identical to Section 
III.c, except we set a uniform density to test the be- 
haviour of frequency as function of time. In Section IILe 
we consider case similar to III.c but with 10 times weaker 
magnetic field such that near the beam injection location 
Lo ce /uj pe (x = 0) = 0.094, unlike in the rest of the paper 
where u ce /uj pe (x = 0) = 0.94 (which is commensurate 
to solar coronal conditions). This is to prove that the 
present study is indeed modelling the situation relevant 
for the type III radio bursts, which emit near the plasma 
frequency, w pe , rather than electron gyro-frequency, ui ce . 



II. THE MODEL 

Wc use EPOCH (Extendible Open PIC Collaboration) 
a multi-dimensional, fully electromagnetic, relativistic 
particle-in-cell code which was developed by EPSRC- 
funded Collaborative Computational Plasma Physics 
(CCPP) consortium of 30 UK researchers. EPOCH uses 
second order accurate FDTD scheme to advance EM 
fields. EPOCH'S particle pusher is based on the Plasma- 
Simulation-Code (PSC) by Hartmut Ruhl, and is a Bird- 
sail and Landon type PIC scheme [29| using Villasenor 
and Buneman current weighting. EPOCH uses a trian- 
gular shape function with the peak of the triangle lo- 
cated at the position of the pseudo-particle and a width 
of twice the spatial grid length. EPOCH utilises the 
Villasenor and Buneman [3(| current calculating scheme 
which solves the additional equation dp/dt — V ■ J to 
calculate the current at each time step. The main advan- 
tage of this scheme is that it conserves charge on the grid 
rather than just globally conserving charge of the parti- 
cles. This means that the solution of Poisson's equation 
is accurate to the machine precision, and when Poisson's 
equation is satisfied for t = it remains satisfied for all 
times. EPOCH has been thoroughly tested and bench- 
marked. 

We use 1.5D version of the EPOCH code which means 
that we have one spatial component along ir-axis and 
there are all three V x ,V y ,V z particle velocity components 
present (for electrons, ions and electron beam). Using 
these, the relativistic equations of motion are solved for 
each individual plasma particle. The code also solves 
Maxwell's equations, with self-consistent currents, us- 



ing the full component set of EM fields E Xl E y , E z and 
B x , By, B z . EPOCH uses un-normalised SI units, how- 
ever, in order for our results to be generic, we use the nor- 
malisation for the graphical presentation of the results as 
follows. Distance and time are normalised to c/u) pe and 
u pe , while electric and magnetic fields to ui pe cm e /e and 
w pe m e /e respectively. Note that when visualising the 
normalised results we use tiq = 10 14 m~ 3 in the dens- 
est parts of the domain, which are located at the left- 
most and rightmost edges of the simulation domain (i.e. 
fix Lo pe — 5.64 x 10 8 Hz radian in the densest regions). 
Here u pe = y / n e e 2 /(eom e ) is the electron plasma fre- 
quency, n a is the number density of species a and all 
other symbols have their usual meaning. We intend to 
consider a single plasma thread (i.e. to use 1.5D ge- 
ometry), therefore space component considered here has 
x = 65000 grid points, with the grid size is Ad/2 for the 
two long runs (Sections III.c and IILe) making maximal 
value for x, x max — 231.170c/w pe , while in the short runs 
(Sections IILa-b and IILd) the grid size is Ad/4 yielding 
the maximal value for x, x rnax = 115.585c/oj pe . Here 
Ad = v t h,e/u pe is the Debye length (v t h,e = y/kT/m e 
is electron thermal speed). Since we would like to re- 
solve full plasma kinetics, our choice of the grid size is 2 
to 4 times better than in Ref.[8| where only spatial grid 
size of IAd was used. Thus, the presented results can 
be regarded as high resolution (sub-Debye length scale), 
guaranteeing a superior capture of kinetic physics. 

We do not fix plasma number density and hence ui pe 
deliberately, because we wish our results to stay gen- 
eral. We demonstrate this on the following example, 
if we set plasma number density to no — 10 14 m~ 3 
(i.e. fix Lu pe = 5.64 x 10 8 Hz radian), this sets Debye 
length at A D = 3.78 x 10~ 3 m=7.11 x 10" 3 c/w pe (us- 
ing T — 3 x 10 5 K). If we set plasma number density 
to no = 10~ 5 m~ 3 , this sets Debye length at Ad = 
1.20 x 10 7 m=7.11 x 10~ 3 c/cl! P£ . Thus, appropriately ad- 
justing plasma number density no, physical domain can 
have arbitrary size e.g. Sun-earth distance (but then un- 
realistically low density has to be assumed). Background 
plasma in our numerical simulation is assumed to be 
Maxwellian, cool T = 3 x 10 5 = const with parabolically 
decreasing density gradient, along the uniform magnetic 
field line. The latter mimics a field line that connects 
Sun to Earth. The only physical parameters that should 
be regarded as fixed are the temperatures of the back- 
ground plasma, that of the electron beam and magnetic 
field (B 0x = 30 gauss). These are set to, plausible for the 
type III bursts values, T = 3 x 10 5 K for the background 
plasma and Tb = 6 x 10 6 K for the beam. This fixes re- 
spective electron thermal speeds to v t h.e = 7.11 x 10~ 3 c 
and Vth,b = 3.18 x 10 _2 c. If an attempt is made to inter- 
pret some type III burst observations, one should keep 
in mind whilst density and hence uj pe can be regarded as 
variable (arbitrary), T, Tb and Bq x are fixed (model spe- 
cific). Ref.fiOJ have shown that parabolic density profile 
n e (r) describes the electron number density to a good 
approximation within few solar radii, Rq. Generally, 
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n e (r) oc r~ 2 plasma number density profile can be well 
understood based on conservation of mass for a spheri- 
cally symmetric constant speed outflow such as Parker's 
solar wind solution. For large radii, r 3> i?0, practically 
all models predict n e (r) oc r~ s scaling with S being close 
to two, e.g. 2.16 in Ref.[32j or 2.19 in Ref.jajjj]. Therefore, 
to a good approximation, we use following density profile 
for the background electrons (and ions): 

n o( x ) = (( x ~ x ma x/2)/(x max /2 + n + )) 2 + n-, (1) 

where Uq(x) is the normalised plasma number density, 
such that for the left and right edges of the simulation 
domain, x = and x = x max , n (0) = n Q (x max ) 
1 . while in the middle, x = x max /2, the parameters 
n+ = (x ma x/2)x (1-yi -n-)/y/l -n- and rt- = 10~ 8 
where chosen such that no(x) drops 10 -8 times compared 
to the edges. This density profile effectively mimics a 
factor of 10 8 density drop from the corona no = 10 14 
m~ 3 to uau — 10 6 m~ 3 at 1 AU. Because numerically 
most precisely implementable boundary conditions are 
the periodic ones, this density profile represents mirror- 
periodic situation when the domain size is effectively dou- 
bled, i.e. at no (a; = 0) = no (a; = x max ) = 1 while 
no(ir = x max /2) = 10~ 8 . This way "useful" or "work- 
ing" part of the simulation domain is < x < x max /2. 
When cases with the beam are considered we set its fol- 
lowing density profile: 

n b (x) = nbe -[(— — /25)/(^/40)] 8 (2) 

which means that the beam is injected at x = x max /25 
and its full width at half maximum (FWHM) is w 
Xmax/20 (see Figure 4(c) dashed curve). 

As, we impose background magnetic field Bq x = 
30 = const gauss along rr-axis, plasma beta in this 
study, based on the above parameters, is set to /3 = 
2(v thii /cf(u pi /u ci ) 2 = n (0)kT/(B 2 /(2u )) = 1.16 x 
10~ 2 atx = 0. It should be noted that the pressure 
balance in the initial conditions is not kept. There are 
two reasons for this: (i) solar wind is not in "pressure 
balance" and it is a continually expanding solar atmo- 
sphere solution; (ii) plasma beta is small therefore it is 
not crucial to keep thermodynamic pressure in balance 
(because its effect on total balance is negligible) and the 
initial background density stays intact throughout the 
simulation time (see e.g. Figure 4(c), thick solid curve). 

EPOCH code allows to set an arbitrary number of 
plasma particle species. Thus, since we intend to study 
spatially localised electron beam injected into the in- 
homogeneous or homogeneous Maxwellian electron-ion 
plasma, we solve for three plasma species electrons, ions 
and the electron beam. The dynamics of the three 
species, which mutually interact via EM interaction, can 
be studied independently in the numerical code. Velocity 
distribution function for electrons and ions is always set 
to 

f e i = e ~(pl+P 2 v+pl)/( 2m ",' kT \ (3) 



where the momenta components, p x ,p y ,p z , include the 
correct electron and ion masses which are different by 
the usual factor of mi/m e = 1836. When cases with the 
beam are considered we set its following distribution 

f b = n b e~ ((Px " Pxo)2+(P! ' _P! ' o)2+p ' )/(2m ' ! '*' cT) . (4) 

where hb is normalised beam number density (fib — 
Ub/rieO ) and it is hb = 10~ 3 throughout this study. This 
choice is on the limit of available current computational 
facilities used - 64 Dual Quad-core Xeon = 64x8 = 512 
processor cores with 4 Tb of RAM. Typical run takes 28 
hours on 512 processors. We have used 6.5 x 10 8 elec- 
trons, 6.5 x 10 s ions and 6.5 x 10 5 beam electrons giving 
a total of 1.30065 x 10 9 particles in the simulation, with 
65000 spatial grid points this means that we sample elec- 
tron and ion phase space very well with 10000 per sim- 
ulation cell. With hb — 10~ 3 this means that globally 
(on average) we only have 10000 x 10~ 3 = 10 electrons 
to represent the beam. Thus, we cannot consider a real- 
istic hb = 10 -5 — 10~ 7 commensurate to type III bursts. 
However as can be seen from Figure 4(c), dashed curve, 
the beam is quite localised, about l/10th of the domain 
length (i.e. twice the FWHM w 2x max /20 = x max /10). 
Thus, in reality 6.5 x 10 5 beam electrons are loaded 
into 65000/10 = 6500 cells providing reasonably good 
6.5 x 10 5 /6500 = 100 beam particles per cell. 

III. RESULTS 

Below we present numerical simulation results for the 
four runs. We use the beam injection initial momentum 
components p x o and p y o to control what type of waves 
can be excited by the beam, as well as study the effect 
of the background plasma density gradient on the gener- 
ation and properties of EM waves. 



A. Inhomogeneous plasma without an electron 
beam 

In this section we present an equilibrium test run where 
the above described initial conditions are evolved for 
§0uj~z without imposing an electron beam. The results 
are presented in Figure 1 where we show time-distance 
plots for the electrostatic (longitudinal to both back- 
ground magnetic field and density gradient) electric field 
E x 1(a) and associated density perturbation, n e — n e o 
1(b); and two components of the transverse electromag- 
netic fields: (E y 1(c), B z 1(d)) and (E z 1(e), B y 1(f)). 
We gather from Figures 1(a) and 1(b) that only low level 
noise (with amplitudes ~ 9 x 10~ 5 for E x and ~ 3 x 10 ~ 2 
for n e — n e o) is generated. This can be attributed to 
so called "shot noise" that is normally present in PIC 
simulations. Figures l(c)-l(f) demonstrate that in the 
electromagnetic emission component also low level (with 
amplitudes few 10 -6 ) drift EM wave noise is generated. 




This can be evidenced by the fact that the slope of bright 
and dark strips is roughly the speed of light. Hence 
the perturbations are travelling with the speed of light. 
Note that the perturbations are generated in all parts of 
the density gradient but they are more prominent in the 



densest parts of the simulation domain, because their am- 
plitudes are also expected to be largest there. No regular 
perturbations of longitudinal magnetic field {B x — Bq x ) 
are found (not shown here). Thus, we conclude that the 
equilibrium without the electron beam is fairly stable 
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(apart from the low level EM drift wave noise). 



B. Inhomogeneous plasma with electron beam 
injected along the magnetic field (9 = 0°) 

In this section we present the results when we inject 
the electron beam along the background magnetic field 
(6 = 0°), where 9 is the beam pitch angle (the angle be- 
tween the initial beam velocity/momentum vector and 
the direction of background magnetic field). Here in 
Eq.(4) at t = we set po x = 0.5jm e c and po y = 0. 
Note that in all cases with the injected beam we solve 
the initial value problem, i.e. electron beam initial drift 
momentum is applied only at t = - we do not re-inject 
the beam at every time step. The results are shown in 
Figure 2 where physical quantities shown are similar to 
that of Figure 1 except for 1(b) where instead we now 
present the time-distance plot for the electron beam, rib- 
The reason why we have chosen to trace the dynamics 
of nb rather than full electron number density perturba- 
tion (background electron population plus the beam) is 
because as it was show in Section 2.3 of Ref.[8], when the 
beam is relatively dense (nb/n e ss few 10~ 2 ), the electron 
number density perturbation is dominated by a wake cre- 
ated by the beam. In Ref . [§] it was shown that when an 
electron beam, with the properties similar to considered 
here, is injected perpendicular to the background mag- 
netic field, the beam excites electrostatic, standing waves, 
oscillating at local plasma frequency, in the beam injec- 
tion spatial location. Here physical situation is differ- 
ent in that now the beam is injected along the magnetic 
field which is more plausible for the type III radio bursts. 
However, surprisingly we still see in Figure 2(a) the simi- 
lar effect, that the standing ES waves are generated in the 
beam injection location. We can estimate the oscillation 
frequency by counting bright yellow strips in the region 
2c/w pe < x < &c/uj pe which is seven starting from the 
first strip. The time elapsed is 50^"^, thus 50/2-7T = 7.96 
and the conclusion is that this standing wave has approx- 
imately the plasma frequency. Note that the small mis- 
match is due to the fact that oj pe in all normalisations 
is taken at the edges of the simulation domain where 
n o(0) = ngf^mai) = 1, while the beam injection location 
is centred on x max /25 where n Q (x max /25) — 0.846. (The 
beam spatial spread here is within 2c/co pe < x < 8c/w pe ). 
In Figure 2(a) 10c/w pe < x < 35c/uj pe we see series of 
oblique strips which is ES wake created by the beam. 
This can be evidenced by the fact that slope of these 
oblique strips is (35 — 10)(c/u! pe )/(50/uj pe ) — 0.5c which 
coincides with the electron beam injection speed. Simi- 
lar conclusion is reached from Figure 2(b) as well where 
the inferred slope of the beam is the same (0.5c). Note 
that for the times t > 30^^/ there is small dip formed 
on the top of the beam (see Figure 2(b)). However the 
beam seems to stay intact which would be expected in 
the quasilinear theory (due to so called beam free stream- 
ing). This is despite the fact that rib/n e — 10 -3 and the 



criterion of weak turbulence regime of quasilinear theory 
e = nbm e vl / (nom e v^ h e ) <C 1 [H[ is actually not met: 
In our case e = 4.94 and also the quasilinear relaxation 
time, r, (time of establishing the plateau in the elec- 
tron longitudinal velocity distribution function) is given 
by r = n e / (nbUJpe) (e.g. [13| ) . In our cause r = lO 3 ^^ 1 . 
Thus, no substantial quasilinear relaxation is expected to 
take place within 50uj pi }. We gather from Figures 2(c)- 
2(f) as in Section III. a that only low level EM drift wave 
noise is generated and we see no generation of regular 
EM waves, commensurate to type III bursts. As dis- 
cussed in the Introduction section when the beam pitch 
angle is zero then no EM waves should be generated; and 
it is only for the oblique pitch angles the EM emission 
generation is possible because then Vb ■ E±_ =/= 0. 



C. Inhomogeneous plasma with electron beam 
injected obliquely to the magnetic field (9 = 45°) 

In this section we present results when we inject the 
electron beam obliquely to the background magnetic field 
(6 = 45°). Here this is achieved by setting po x — 0.5jm e c 
andpoy = 0.57m e c in Eq.(4) at t = 0. Note that since we 
intended to consider a numerical run for a twice longer 
time (lOOWyg 1 ) than in other above Sections, we have dou- 
bled the spatial domain size whilst keeping the same to- 
tal number of particles 1.30065 x 10 9 . This implies that 
now the spatial grid size is Ad/2, not Ad/4 as in other 
above Sections. We gather from Figure 3(a) that again 
in the beam injection spatial location standing ES wave 
oscillating at local plasma frequency is excited. However, 
the ES wake of the beam, oblique yellow lines between 
30c/w pe < x < 65c/w pe , detaches from the standing ES 
wave and becomes localised. The ES wake also travels 
with the correct speed of ~ 0.5c. Similar conclusions 
can be reached from analysing Figure 4(a) which reveals 
more detailed spatial structure of the ES wave oscilla- 
tion and the ES beam wake. Figure 3(b) corroborates 
the beam travel speed of 0.5c as well as reveals minor 
deviation from the quasilinear theory free streaming, by 
appearance of spikes on the top of the beam. Figures 
3(c)-3(f) present time-distance plots for the two compo- 
nents of the transverse electromagnetic fields: (E y 3(c), 
B z 3(d)) and (E z 3(e), B y 3(f)). We learn that as the 
beam pitch angle now is 6 = 45° ^ escaping EM ra- 
diation is generated. In the beam injection spatial loca- 
tion, 4c/w pe < x < 15c/w pe , we see strong interference 
pattern between the standing (trapped) ES and escap- 
ing EM radiation. The EM waves travel in both direc- 
tions, and because of the periodic boundary conditions, 
waves that travel to the left, appear on the right side 
of the simulation domain (x > 150c/o;p e ). Figure 4(c) 
shows normalised electron number density at t = lOOcjpg 1 
(thick solid curve), beam spatial profiles at t — (dashed 
curve) and t = lOOo;^ 1 (thin sold curve). Note that 
when plotting the spatial profiles of nb, we scale it by 
a factor of 10 3 so that it is clearly visible. We gather 
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that the background electron population number density 
stays unchanged throughout the simulation, compared 
to t = 0. This serves as additional proof that our initial 
conditions without the beam are stable. Also, compar- 



ing Figures 4(a) and 4(c) we confirm that indeed the 
beam and ES wake travel the same distance at the same 
speed of « (65 — 15)c/ujp e /(l00/ujpe) = 0.5c. Figures 
4(b) and 4(d) show a more detailed spatial structure of 



E, (a) 



n b (b) 







1 

50 


100 

* [c/Up] 


150 


200 







100 150 

» [e/»(J 



2.OB-03 O.Oe-l-OO 



m 



9.9a-D+ 2.08-03 3.08-03 

B 2 (d) 



I 



50 100 150 200 



+.28-04 -3,58-04 



-a.ae-05 1.58-04 

B„ (f) 



w 



I 



1 



1 , 




, I 1 , 

50 


r 

100 

H [c/UpJ 


150 


200 





50 


100 


150 


200 





+.1B-04 -3,28-04 



FIG. 3: As in Figure 2 but for the beam pitch angle 9 = 45°. This figure pertains to Section III.c. 



the transverse, escaping EM radiation components. We 
gather that these actually consist of two parts: (i) the 
part within < x < 65c/aj pe corresponds to the non- 
escaping ES standing waves and the ES wake of the beam; 
and (ii) the part within 65c/cj pe < x < 115c/w pe that 



corresponds to the escaping EM radiation. Figures 4(e) 
and 4(f) show electron (including the electron beam) and 
ion longitudinal velocity distribution function time evolu- 
tion. The considered momentum range is 2iVth,i^ni, which 
is then converted to the velocities (in the relevant Fig- 
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ures) by using 3vth,inii / '(7 TO e,ic) . We gather from Figure 
4(e) that background electron population distribution re- 
mains unchanged (dashed it = 0) and solid (t — lOOo;^ 1 ) 
curves, centred on v x — do overlap to a plotting accu- 
racy), while the electron beam starts to show a tendency 
of plateau formation, according to the quasi-linear the- 
ory. Note that in this Section simulation end time is 
100 uj~ l } while the quasi-linear relaxation time is 1000 
u!p L } . Recall however that we are not strictly speaking in 
the quasi-linear regime because e « 5. The ion distri- 
bution function (Figure 4(f)) shows no noticeable by-eye 
change. Thus despite the fact that ions in the simula- 
tion are treated as mobile, the ion population shows no 
dynamics in the velocity space. 

Next we attempt to produce a synthetic (simulated) 
dynamical spectrum. We take a snapshot of the spa- 
tial profile of one of the transverse EM components, 
E y (x,t = lOOojpg 1 ), and cast it into the temporal 
dependence by putting x = ct. Thus in Figure 5(a) we 
see the same pattern in E y as in Figure 4(b) but now it 
appears as function of t normalised to oj^ e . Note that 
we do not include EM wave which appear on the right 
due to the periodic boundary condition, i.e. we restrict 
ourselves to the range x = [0, 128c/cj pe ] (i.e. the same 
as t = [0,128^"/]). We then generate a wavelet power 
spectrum for the E y (t). Wavelet software was provided 
by C. Torrence and G. Compo, and is available at URL: 
http : //atoc . Colorado . edu/research/wavelets/. 
We gather from Figure 5(b) that in the time interval 
< x < 65w pe 1 the wavelet power spectrum is fiat 
(period/frequency does not change in time). This corre- 
sponds to ES oscillation part, oscillating at the plasma 
frequency (with a prefactor of 2tt). In the time interval 
65 < x < 128^"/ wavelet power spectrum corresponds 
to escaping EM radiation part and we clearly see a de- 
crease of EM signal frequency in time. One has to realise 
however that in the Figure 5(b) it appears that large 
period Plf = 15/(27r) (low frequency) at t 65lj~J- 
shows up first and then low period Phf = 9/(27r) (high 
frequency) at t w 115^"/ follows. This has a simple 
explanation that we obtained the time series of E y by 
putting x — ct. In reality for a distant observer located in 
a point, the high frequency (high density) would appear 
first followed by a low frequency radio signal. Note also 
that the frequency decrease nicely follows the number 
density decrease. In Figure 4(c) dashed line peak which 
represents the beam at t — is located at the normalised 
background electron number density (thick solid line) 
of 0.846. By the time beam reaches its final destination 
by time t — lOOo;^ 1 (thin solid line) the density has 
dropped to ~ 0.3. Therefore we would have expected 
that the corresponding plasma frequency has dropped 
by a factor ^/0. 3/0. 846 = 0.595. Indeed we gather from 
Figure 5(b) that Phf /Plf = 9/15 = 0.6. In other 
words, ujpxp /^Jp,HF = \J n LF / n H f = ^0.3/0.846 = 
0.595 w Phf/Plf = 9/15 = 0.6, where the notation 
is straightforward. Therefore we conclude that the 
frequency decrease in the synthetic dynamical spectrum 



is commensurate to the plasma frequency (electron 
number density) decrease along the beam propagation 
path. Note that it is only in the interplanetary type 
III radio bursts the frequency drops by many orders of 
magnitude. However, the bursts that occur in the solar 
corona show frequency drops by about a factor of two 
(as in our simulation), are not uncommon. For example 
decameter type II bursts which have a fine structure 
in the form of type Ill-like bursts. The drift rates of 
these sub-bursts are close to the ordinary type III bursts 
velocity, but their duration is essentially lesser [34| . 



D. Homogeneous plasma with electron beam 
injected obliquely to the magnetic field (6 = 45°) 

In this section we present the results when we inject 
the electron beam obliquely to the background magnetic 
field (9 — 45°) as in Section III.c. However, contrary 
to the Section III.c here we consider a plasma with uni- 
form background number density uq = 1. We see from 
Figure 6(a) that as in Section III.c in the location where 
the beam was injected standing ES waves are generated, 
oscillating at local plasma frequency. However, as better 
seen from Figure 7(a) now the ES wake created by the 
beam does not have enough time to detach itself from 
the standing wave. This is due to the fact that this is a 
shorter run (t enc t — 500^). We also gather from Figures 
6(b) and 7(c) that the beam travels the correct distance, 
commensurate to its speed. Figures 6(c)-6(f) show time 
distance plots of transverse EM components that are gen- 
erated by the Langmuir (ES) waves. Their more detailed 
spatial profiles are shown in Figures 7(b) and 7(d). It is 
evident by eye that in the uniform plasma case there is 
no frequency decrease with time for the generated EM 
components. Therefore we do not produce the synthetic 
dynamical spectrum as in Section III.c. Figures 7(e) con- 
firms that the predictions of the quasilinear theory, that 
by time t — 500;"^ we see (i) no noticeable quasi-linear 
relaxation because the plateau is expected to develop in 
quasi-linear relaxation time of 1000 U!~J-; (ii) electron free 
streaming is also evident. No noticeable change in ion 
velocity distribution function is seen either. Note that 
the escaping EM radiation is generated in the uniform 
plasma density case too. This indicates that the density 
gradient plays no role in the EM emission generation. 



E. Inhomogeneous plasma with electron beam 
injected obliquely to the magnetic field (6 = 45°), 
weak magnetic field case 

In this Section we consider the case similar to III.c 
but with ten times weaker magnetic field, B = 3 gauss, 
such that oj ce /ujp e (x = 0) = 0.094, unlike in the rest of 
the paper where uj ce /u> pe (x = 0) = 0.94 (which is more 
appropriate to solar coronal conditions). This is to un- 
ambiguously demonstrate that the present study is in- 
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deed relevant for the type III radio burst emission, via 
the 1.5D non-zero pitch angle (non-gyrotropic) electron 
beam quasilinear relaxation, and subsequent emission 
at the plasma frequency rather than the electron gyro- 
frequency emission. In the laboratory plasma there are 
microwave generation devices, such as Gyrotron, in which 
EM radiation is produced at electron gyro-frequency, 
uj ce . At first sight, it would seem probable that since for 
the considered model parameters throughout this paper 
uj ce /uj pe {x = 0) = 0.94, what we report is the Gyrotron 
type EM radiation. The issue can be settled by e.g. low- 
ering the magnetic field value. The results are presented 
in Figs. 8-10, which are mirror analogs to Figs. 3-5, ex- 
cept that here uj ce /uj pe (x = 0) = 0.094. We gather from 
Fig. 8a that the generated ES (Langmuir) component be- 
haviour is nearly identical to that of Fig. 3a, i.e. we again 
see the generation of standing ES waves, oscillating at 
plasma frequency, u) pe , in the beam injection spatial lo- 



cation, 4c/ujp e < x < 15c/w pe . The electron beam dy- 
namics is also identical (cf. Fig. 8b and Fig. 3b). There 
are notable differences in the generated transverse to the 
magnetic field EM components: By comparing Figs. 8c- 
8f to Figs.3c-3f, we see that in the escaping EM ration 
there is no interference pattern, i.e. there is no inter- 
ference between the standing (trapped) ES and escaping 
EM radiation. This is because now ES oscillation no 
longer appears in E y , E z , B y , and B z transverse EM 
components. What is crucial that in the time-distance 
plots, Figs.8c-8f, we still see the same number of bright 
(and dark) strips, which is indicative of the fact that the 
escaping EM radiation oscillates again at approximately 
plasma frequency, to pe , and not at electron cyclotron fre- 
quency UJ ce . 

We gather from Fig. 9, which corresponds closely to 
Fig. 4, but here uj ce /uj pe {x = 0) = 0.094, that most of the 
conclusions reached when considering Fig. 4 apply here 
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FIG. 5: (a) Time evolution of E y (see text how spatial x- 
dependence was cast into time f-dependence), (b) wavelet 
power spectrum of (a). This figure pertains to Section III.c. 



too. However, a notable difference is that by comparing 
Figs. 4b and 9b we deduce that in transverse electric field 
component E y , for < x < 65c/w pe , we no longer see the 
standing (trapped) ES component, i.e. only the escaping 
EM radiation is present. 

Fig. 10 has been produced in the same way as Fig. 5, 



except that now uj ce /u> p 



0) = 0.094. We gather 



from Fig. 10 that, contrary to Fig. 5 in which both the 
ES oscillation and escaping EM component were present, 
now we see only escaping EM wave which clearly shows 
the drift towards lower frequencies. We also observe that 
as in Fig. 5, oj p .lf/^p,hf ~ 0.6 that is commensurate to 
the background plasma density decrease. 



IV. CONCLUSIONS 

A quote from [35], p. 95, summarises the state of the 
matters in theoretical understanding of solar type III ra- 
dio bursts rather well: "Our understanding of plasma 
emission is in an unsatisfactory state. It seems that the 
problems with our understanding of plasma emission are 
of an astrophysical nature and will eventually be solved 
through new observational data. There are several differ- 
ent possible mechanisms which can lead to fundamental 
plasma emission and it is still not clear which is the rel- 
evant one in practice. Although this leaves the theory of 
fundamental plasma emission in a somewhat uncertain 
state, the theory for second harmonic is well understood; 
there seems to be no reasonable alternative for the coales- 
cenceprocess L + L — > T" . Also Introduction section of 
Ref. 9] provides a good, critical overview of the possible 
mechanisms which generate the type III burst EM radi- 
ation. These include, (i) the classical plasma emission 
mechanism that is based on non-linear wave-wave inter- 
action between Langmuir, ion-acoustic and EM waves; 



(ii) a linear mode conversion, in which almost monochro- 
matic Langmuir z-mode interacts with the density gra- 
dient, partly reflecting and partly converting into the 
EM radiation; (iii) the quasimode mechanism in which 
forward- and backward- propagating Langmuir waves 
generate a quasinormal electrostatic mode at 2u> pe which 
further converts into EM harmonic radiation and (iv) the 
antenna radiation model which involves direct radiation 
of charged particles that oscillate at uj pe and drive cur- 
rents at 2uj pe . In this work we presented the results which 
show that 1.5D non-zero pitch angle (non-gyrotropic) 
electron beam can also produce escaping EM radiation 
at w pe which seems to successfully mimic the observed 
solar type III radio bursts. A clear distinction needs to 
be drawn between the best studied type III burst mech- 
anism, the classical plasma emission [36] , and the model 
presented in this work. In the plasma emission mech- 
anism non-linear wave-wave interaction between Lang- 
muir, ion-acoustic and EM waves requires that the beat 
conditions u>i + uji = lj and k\ + lt2 — k to be satisfied. 
The emission formula (i.e. the three wave interaction 
probability) for the relevant process L + s — > T (coales- 
cence of Langmuir and ion-sound wave to produce trans- 
verse EM wave) includes a cross vector product factor 
\k L x k T \ 2 (see e.g. Eqs. (26.24) and (26.25) from Ref. (37|. 
This implies that, whilst electron beam and Langmuir 
turbulence dynamics can be treated in ID spatial dimen- 
sions (and there is large body of work that deals with 
the ID quasilinear theory), the correct treatment of act 
of EM emission needs 2D spatial dimensions. This is 
because in ID case the factor \k L x k T \ 2 is identically 
zero, because the angle between wave vectors of Lang- 
muir wave, kh, and EM wave, kr, is zero. This has to be 
distinguished from the pitch angle 9 which in our notation 
is the angle between the particle (electron) beam injec- 
tion direction and the background magnetic field. We re- 
mark however that the plasma emission mechanism equa- 
tions use "random phase approximation" (see e.g. p. 383 
from Ref. [37]). This is because the extraneous, quadratic 
non- linear current (see their Eq.(26.2)) depends on the 
phases of the beating fields and some assumption needs 
to be made concerning the phase. A priori, it is not 
at all clear that in the case of our situation, in which 
non-zero pitch angle (non-gyrotropic) electron beam is 
injected, the phases are random. Hence whether the fac- 
tor \kL x kx\ 2 is applicable in our case. Without further 
in-depth analysis it would be safe to conclude that our 
simulations do not involve the classical plasma emission 
processes. 

We have performed-high resolution (sub-Debye length 
grid size and 10000 particle species per cell), 1.5D 
Particle-in-Cell, relativistic, fully electromagnetic simu- 
lations to model electromagnetic wave emission genera- 
tion in the context of solar type III radio bursts. We 
studied the generation of EM waves by injecting a super- 
thermal, hot beam of electrons into a plasma thread 
that contains uniform longitudinal magnetic field and a 
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to Section Ill.d. 



parabolic density gradient along the magnetic field. We 
have considered five cases: 

(i) As an initial equilibrium test, we find that the phys- 
ical system without electron beam is stable and only low 



amplitude level electromagnetic drift waves (noise) are 
excited. 

(ii) The beam injection direction is then controlled by 
setting either longitudinal or oblique initial electron drift 
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speed/momentum, i.e. we set different beam pitch an- 
gles. In the case of zero beam pitch angle, i.e. when 
Vb ■ E± = 0, the beam excites only ES standing waves, 
oscillating at local plasma frequency. This oscillation oc- 
curs strictly in the beam injection spatial location and 
only low level electromagnetic drift wave noise is present 
(no regular EM waves are generated by the beam). 

(hi) In the case of oblique beam pitch angle, i.e. when 
Vb ■ Ej_ ^ again ES waves with similar properties arc 
excited. In this case however because the beam can in- 
teract with the EM waves, it generates EM waves with 
the properties commensurate to type III radio bursts. 
In particular, wavelet analysis of transverse electric field 
component shows that as the beam moves to the regions 
of lower density and hence lower plasma frequency, EM 
wave frequency drops accordingly. 

(iv) When we remove the density gradient, an electron 
beam with an oblique pitch angle still generates the EM 



radiation, but now no frequency decrease is produced. 

(v) In order to prove that the generated, by the non- 
zero pitch angle beam, EM emission oscillates at the 
plasma frequency, we also consider a case when the mag- 
netic field (and hence the cyclotron frequency) is ten 
times smaller. 

Using fully kinetic plasma model our results also 
broadly confirm (i) the fact that in order to excite es- 
caping EM waves, the electron beam should have non- 
zero pitch angle; i.e. there should be a non-zero pro- 
jection of the electron beam injection velocity vector on 
the transverse EM electric field vector, (ii) quasilinear 
theory predictions, namely quasilinear relaxation time- 
scale and free streaming assumptions were corroborated 
via fully kinetic simulation, in a realistic to the type III 
burst magnetic field geometry, (iii) The observational 
fact that there should be a EM emission frequency drift 
in time in the inhomogeneous plasma case has been also 
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FIG. 8: As in Figure 3 but for 10 times weaker magnetic field such that uj ce /uj pe (x = 0) = 0.094. This figure pertains to Section 
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confirmed via production of the simulated (synthetic) dy- 
namical spectrum for the first time. 

The presented model will be used in the future for 
the forward modelling of the observed type III burst dy- 



namical spectra (2D radio emission intensity plots where 
frequency is on y-axis and time on ir-axis). The main for- 
ward modelling goal will be inferring the electron number 
density profile along the beam propagation paths. 
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FIG. 9: As in Fig.4 but for 10 times weaker magnetic field such that ui C e./upe.(x = 0) = 0.094. This figure pertains to Section 
IILe. 



We would like to close by pointing out some pertinent 
limitations of the considered model in its direct applica- 
bility to the solar type III radio burst observations. The 
issues are: 

(i) In the solar type III bursts, electron beams propa- 
gate large distances without being depleted by the gen- 
erated Langmuir waves due to bump-on-tail instabil- 
ity (also called beam-plasma instability). In the lin- 
ear regime, the timescale for the quasilinear relaxation, 
t = n e / (ribUjpe) and hence timescale of the beam deple- 
tion is prescribed by the ratio of the electron beam and 
background plasma densities. To be precise, r is the in- 
verse of the linear bump-on-tail instability growth rate 
7 = ujp e (nb/n e )(v^/Sv^), where Svb is the electron beam 
velocity thermal spread. Thus for Vb ~ <5i>b, r ~ 1/7. 
Due to the above described computational limitations, 
at present, it was impractical to set rib/n e to the ob- 
served values 1CT 5 — 10~ 8 ; This may affect the pro- 
cess of beam (re-)generation by the time-of-flight effects; 
Also, it is known that electron beam may be stabilised 
by non-linear effects [38j |. In the non-linear stimulated 
scattering processes, the wavenumbers are drawn out of 



the resonance. This leads to energy transfer rate be- 
tween the beam electrons and Langmuir wave at a much 
slower rate than quasilinear relaxation time, that effec- 
tively leads to the non-linear stabilisation of the bump- 
on-tail instability. For the parameters commensurate to 
the type III bursts, the condition for the stabilisation, 
rib/n e <C v the/ v b ^ (n r ie/fni)(dvb/vb) is met in most 
cases. Therefore the non-linear stabilisation is likely to 
play a major role. (e.g. Ref.[39j], pp.184-187 or [Hj]) 

(ii) The spatial scale of the density gradient (decrease 
of ujp e by a factor of 10 4 over 65,000 Debye length) is 
not realistic (we have to remember that in PIC sim- 
ulations one uses a smaller than in reality number of 
"super-particles" and not real electrons and protons). 
This may affect criteria for onset of instabilities caused by 
plasma density inhomogeneities; For example, when the 
electron beam moves along the density gradient, Lang- 
muir wave phase velocity will change while the beam ve- 
locity remains constant when there is no strong relax- 
ation [40]. If the instability growth rate is much less 
than the reciprocal of the time of escape from the res- 
onance, the beam stabilises as it no longer loses energy 
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FIG. 10: As in Fig. 5 but for 10 times weaker magnetic field 
such that u} ce /u) pe (x = 0) = 0.094. This figure pertains to 
Section IlI.e. 

to the wave generation. The condition for the stabil- 
isation is L < (n e /nb)(3Vt h e /(vbU) pe )), where L is the 
characteristic spatial scale of the plasma density inho- 
mogeneity (see e.g. Ref . [4ll] . p. 119). For the broad 
range of the solar coronal conditions as well as for the 
set of parameters considered in this paper (from Figure 
4c, thick sold curve, we see that plasma density drops by 
a factor of 2 over a length scale of L ps 30c/u> pe , whereas 



(n e /nb){3vt h e /(vbUj pe )) — 0.3c/oj pe ) the latter inequality 
is not met. 

(iii) Moreover the observed electron beam pitch angles 
are also much smaller [42| than considered in the present 
model. Despite these limitations, the present model pro- 
vides a proof-of-concept for the EM emission generation 
in the context of type III solar radio bursts. 

On a positive note, the considered regime may pro- 
vide an important diagnostic to laboratory laser plasma 
or thermonuclear fusion studies as in both cases non- 
thermal beams of electrons are frequently present. 
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